Debye-Hiickel theory for spin ice at low temperature 
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At low temperatures, spin ice is populated by a finite density of magnetic monopoles — pointlike topologi- 
cal defects with a mutual magnetic Coulomb interaction. We discuss the properties of the resulting magnetic 
Coulomb liquid in the framework of Debye Huckel theory, for which we provide a detailed context-specific 
account. We discuss both thermodynamical and dynamical signatures, and compare Debye Huckel theory to 
experiment as well as numerics, including data for specific heat and AC susceptibility. We also evaluate the 
entropic Coulomb interaction which is present in addition to the magnetic one and show that it is quantitatively 
unimportant in the current compounds. Finally, we address the role of bound monopole anti-monopole pairs 
and derive an expression for the monopole mobility. 
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I. INTRODUCTION 

Spin systems with long-range interactions, where each spin 
interacts with all others, present a formidable challenge to the- 
oretical analysis. While simplifications occur in the limit of 
infinite range interactions, the case of dipolar interactions in 
three spatial dimensions is particularly complex due to their 
(non-integrable) algebraic decay combined with angular de- 
pendence on the spin direction'. As the determination of the 
behaviour of even a spin model with only short ranged com- 
peting interactions can pose a non-trivial problem, it is a priori 
not obvious how long-range interactions can be treated. 

A remarkable counterexample to this case for pessimism 
is provided by spin ice^, a dipolar Ising magnet on the py- 
rochlore lattice that fails to order down to the lowest temper- 
atures accessed. To a fine approximation, which we detail 
below, spin ice is governed by a model dipolar Hamiltonian 
about which quite a lot is known. 
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where J™ is the exchange interaction truncated at the nearest- 
neighbour level, the spins Si point parallel to the local [111] 
axis (see Fig. [T]), and /iq is the vacuum permeability. The 
rare earth spins S^ have typically a dipole moment of approx- 
imately 10 i-iB {jJ-B = Bohr magneton). 

Most prominently, the model Hamiltonian has an extensive 
set of ground states which can be specified by a purely local 
"ice rule". Their entropy is known to an excellent approxima- 
tion due to Pauling's work already in the context of water ice 
and it has been observed experimentally-'. The T — ^ static 
coiTelations are averages over this ground state manifold and 
their long distance forms are known as they are described by 
an emergent gauge field in the Coulomb phase^**, which have 
also been observed experimentally^ " . 

At low temperatures the physics of the system turns out 
to allow a further simplification. The excitations about the 
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FIG. 1 . The magnetic moments in spin ice reside on the sites of the 
pyrochlore lattice, which consists of corner sharing tetrahedra. These 
sites are at the same time the midpoints of the bonds of the diamond 
lattice (black) defined by the centres of the tetrahedra. The Ising 
axes are the local [111] directions, which point along the respective 
diamond lattice bonds. The bonds of the pyrochlore lattice are in the 
[110] directions, while a line joining the two midpoints of opposite 
bonds on the same tetrahedron defines a [100] direction. 



ground state manifold take the form of magnetic monopoles — 
pointlike defects that interact via a magnetic Coulomb inter- 
action energy which is independent of the background spin 
state'^. In this regime, the magnetic monopoles are sparse, 
as their number is suppressed on account of their excitation 
gap. This in turn has two implications. Firstly, the static 
correlators continue to be dominated by their known T — 
forms up to the inter-monopole separation, whereupon they 
match onto the asymptotics of the paramagnetic phase- . Sec- 
ondly, the low temperature thermodynamics of spin ice can 
be transformed from that of a dense set of localised dipolar 
spins to that of a dilute set of itinerant Coulombically inter- 
acting particles — a (magnetic) Coulomb liquid as first noted 



inRef.lil 



H 



47r ^ ry 



AE 



2fi/ad 



where the charges qi take the values ±2/i/ac(, M — IO/jb be- 
ing the dipole moment of a spin and a^ the distance between 
the centres of adjacent tetrahedra (diamond lattice constant in 
Fig.[T]l, and A is the energy cost of a monopole. 

The transformation is extremely helpful as much is known 
about Coulomb liquids, with a venerable history spanning 
fields from statistical physics all the way to the chemistry of 
electrolytes. Indeed, the known properties of the Coulomb 
liquid have led to an explanation of the iiquid-solid' phase 
transition of spin ice in a [111] fieldJ^, as well as of its mag- 
netic specific heat'" in zero field. More recently, much atten- 
tion has been devoted to the study of the "magnetricity"'"* in 
these "magnetolytes"'^, the equilibrium and non-equilibrium 
behaviour of such a magnetic Coulomb liquid, inspired by the 
analogous electric phenomena such as the Wien effecti^ii^. 

In this paper, expanding on our previous work in Ref. [l^, 
we develop a low-energy theory for spin ice in the framework 
of the Debye-Hiickel (DH) theory of a dilute Coulomb liq- 
uid. DH theory will be familiar to readers from many different 
discipUnes but to our knowledge has never been applied to a 
three-dimensional magnetic material before the advent of spin 
ice. 

The purpose of this paper is two-fold. First, it gives a de- 
tailed and context-specific account of the DH theory for spin 
ice. Second, its ability to model experimental data is un- 
derlined. In particular, we show that an existing framework 
to describe the dynamics of spin ice, when supplemented by 
DH theory, provides improved agreement with existing exper- 
imental and numerical data on the AC-susceptibility of spin 

This is perhaps as good a point as any to digress and address 
the concerns of readers who may be worried that our replace- 
ment of spins by monopoles is too good to be true. Here three 
points are in order. First, as we have already noted above, the 
spins do enter the static correlations but in a manner that is 
understood. Second, a given monopole configuration can be 
"dressed" by many spin configurations. However summing 
over these dressings generates an effective entropic Coulomb 
attraction between the monopoles at long wavelengths (see 
e.g., Ref.l22l) which can also be included in the Coulomb/DH 
framework. We will address this point is Sec. |V] and find that 
the entropic effect can be ignored for the present set of spin 
ice compounds. Third, there is still a remaining issue that 
not all monopole configurations are in fact compatible with 
some spin configuration, and moreover the spins can induce 
non-trivial structure to the monopole energy landscape which 
in turn can significantly alter dynamical properties of spin ice 
out of equilibrium^''. However, these are weak constraints on 
the Coulomb framework and it seems highly unlikely that they 
play any role in determining equilibrium properties. 

We close the introduction by remarking on the range of ap- 
plicability of the Coulomb liquid/DH theory framework in the 
actual compounds (see Fig.|2]i. At high temperatures, above a 



scale Tp, we are in a conventional paramagnetic regime where 
the monopoles are dense. Below Tp the monopoles become 
sufficiently dilute that they can be treated by DH theory. At 
a much lower temperature Td, the Coulomb phase is unsta- 
ble to ordering transitions^"'"^^, the details of which are not 
entirely settled. For the model Hamiltonian, Td = 0. While 
the Coulomb liquid framework should thus apply in the range 
Td < T < Tp, the equilibrium DH treatment runs into prob- 
lems around a temperature Tj > Td where the system falls 
out of equilibrium before any ordering is visible. Much of the 
interest in the spin ice compounds Dy2Ti207 and Ho2Ti207 
derives from the fact that Td,Tf < Tp, so that there is a win- 
dow where Coulomb physics is well visible. 



ordered phase 



spin ice (Coulomb) phase 




out of equilibrium (expm.) 



FIG. 2. Schematic illustration of the different temperature regimes 
in spin ice, separated by Td, Tf, and Tp as explained in the text. 
The putative ordering below Td appears to be prevented by freezing 
of the magnetic degrees of freedom below Tf, as evidenced e.g., by 
a discrepancy between field-cooled and zero-field-cooled magnetisa- 
tion. At temperatures of about Tp, the materials cross over to a trivial 
paramagnetic behaviour. 

The remainder of this paper is organised as follows: we 
first provide DH background, discuss specificities of its ap- 
plication in the spin ice setting, discuss its range of validity 
and finally apply it to experiment. In addition, we discuss 
two other topics of import in this context. Firstly, we deter- 
mine the size of the entropic Coulomb interaction between 
monopoles. Secondly, we compute the low-temperature mo- 
bility of magnetic monopoles in spin ice with a single-spin flip 
dynamics believed to be appropriate for experimental com- 
pounds Dy2Ti207 and Ho2Ti207. 



II. DEBYE-HUCKEL FREE ENERGY 

We now turn to the application of DH theory to spin ice. 
The reader not interested in details of the formalism can skip 
ahead to Section lTV] 



A. Non-interacting monopoles 

To lay the foundation, let us start by considering the sim- 
ple case of non-interacting monopoles, corresponding to a 
nearest-neighbour spin ice model. Since the monopole de- 
scription of spin ice is valid only when the density of defective 
tetrahedra is sufficiently small, i.e., at low temperatures, we 
consider only the less costly defects (3in-lout and 3out-lin 



tetrahedra) and neglect charge 2 excitations altogether (4in- 
Oout and 4out-0in) as they cost four times as much energy. 
The internal energy U of the system is thus proportional to 
the number of monopoles N, 
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(2.1) 



where A is the energy cost of an isolated monopole (assumed 
in the following to be measured in Kelvin) and p = N/Nt is 
the monopole density per tetrahedron. 

The number of configurations that an ensemble of N/2 pos- 
itive (hard-core) monopoles and N/2 negative ones can take 
on a lattice of Nt sites (Nt being the total number of tetrahedra 
in the system) is given by 
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Using Stirling's approximation in the large Nt and large N 
limit, we obtain the S = kslnW 'entropy of mixing', 

5* = S/kB 

= -Nt [2(p/2) In (p/2) + (1 - p) ln(l - p)] (2.3) 

with a concomitant free energy per spin 
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where the number of spins is twice the number of tetrahe- 
dra, Ns = 2Nt. Minimizing with respect to p, we obtain the 
known expression for the total monopole density 
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2exp(-A/T) 
l + 2exp(-A/r)' 



(2.5) 



For small T, and hence small p„n, Pnn — 2 exp(— A/T). For 
large T, Eq. ( I2.5l l tends asymptotically to the value 2/3, which 
is clearly incorrect - as expected since random Ising spins on 
a pyrochlore lattice yield a density Piandom = 5/8 of defective 
tetrahedra. This can be seen e.g., if we consider a single tetra- 
hedron: out of the 2"^ — 16 allowed Ising configurations, only 
6 satisfy the 2in-2out condition and the remaining 10 config- 
urations violate charge neutrality. 



B. Debye-Hiickel contribution 



One of the major approximations in Sec. IIIAl is the fact that 
the long range Coulomb interactions between the monopoles 
were entirely neglected'-. Taking advantage of the analogy 
between spin ice defects and a two-component Coulomb liq- 
uid (in the absence of appplied magnetic fields), we can use 
the Debye approximation to estimate the magnetostatic con- 
tribution to the free energy (in degrees Kelvin per spin)i2i 
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where py = N/V is the dimensionful volume density of 
monopoles and ad is the distance between the centres of two 
neighbouring tetrahedra (i.e., the dual diamond lattice con- 
stant). 

It is convenient to express the dimensionless quantity a^K 
in terms of the Coulomb energy between two neighbouring 
monopoles Enn = /Uo'?^/(47rad fcs). 



QdK = v47r 



Em 



■iPvad) 



(2.7) 



Here q stands for the magnitude of the monopole charge {q = 
2p/ad, where p is the rare earth magnetic moment'-). 

There are 8 diamond lattice sites in a 16-spin cubic unit cell 
of side (4/\/3) Ud- The total volume of the system can then 
be written asV = {Nt/8){4:/V3f a^ and 
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As a result, we arrive at 
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(2.9) 
(2.10) 



where the last equation defines the function a{T). In the 
low temperature limit, the magnetostatic contribution scales 
as p'^/^, namely 
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We can then combine Eqs. (|2.9l l and ( |2.10t with Eq. ( |2.4| i 
from Sec. Ill Al to obtain a mean field free energy - per spin in 
degrees Kelvin - of an ensemble of N monopoles on a lattice 
with long range Coulomb interactions: 
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Note that this reduces to the non-interacting limit if we set 

F„n - 0. 

Minimizing with respect to the defect density p, one obtains 
a self-consistent set of equations: 
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Unfortunately, Eq. ( 12.131 1 cannot be solved analytically and 
one has to resort to numerical methods to obtain p{T). We 
find that the recursive approach 
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converges with acceptable accuracy in less than 5 iterations. 
Substituting p = pi^oo — P5 into Eq. ( 12.12b we obtain nu- 
merically the approximate free energy of dipolar spin ice as a 
function of temperature. 

Between Eqns. ( I2.12l i and ( 12.131 1 we have obtained the free 
energy for monopoles in the DH approximation. From this 
one can compute several thermodynamic quantities of interest 
(see e.g., Sec. lVIAt . 



III. SPIN ICE PARAMETERS AND DH INTERNAL 
CONSISTENCY 

We first derive the parameters describing the Dy2Ti207 and 
Ho2Ti207 spin ices within the dumbbell model'^ in the sub- 
sequent subsection. Following the determination of the pa- 
rameters, we discuss the range of temperatures over which the 
treatment is valid. 



4.6 pLB/k ~ 4.2810-13 J/(Teslam) (see Ref. H^ 
plementary Information therein). 

Inserting the dipolar coupling constant 



and Sup- 



D = 



Pa P- 
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(/io/47r = 10-^ N/A2, kg = 1.38 lO'^^ j/j^) j^^^j t^e nearest- 
neighbour exchange coupling J ~ —3.72 K for Dy2Ti207 
(J ~ —1.56 K for Ho2Ti207) into the expression for the bare 
cost of a single isolated monopole in Ref. |T2, we obtain 



1 2 2J 8 

2 °^ 3 3 



1 



D 



(3.1) 



_ ( 4.35 K for Dy2Ti207 (J = -3.72 K) 
~ 1 5.79 K forHo2Ti207 (J = -1.56K) ' 

The energy of two monopoles at nearest neighbour distance 
is: 
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(3.2) 



Therefore, the creation of two neighbouring monopoles by a 
single spin flip event in a spin ice configuration where all tetra- 
hedra satisfy the 2in-2out rules incurs an energy cost 



A, = 2A - E,, 



5.64 K forDy2Ti207 
8.52 K forHo2Ti207 



(3.3) 



A. Spin ice parameters in the dumbbell model 

The usefulness of the dumbbell model lies in the fact that it 
correctly captures the long-distance form of the dipolar inter- 
action - as well as the magnetic Coulomb interaction between 
the monopoles - while preserving the degeneracy of the spin 
ice states. At the same time, a model of such simplicity can- 
not do justice to the full short-distance structure of the inter- 
actions present in the real compound, which include further- 
neighbour superexchange as well as quadrupolar interaction 
terms between the spins. We will thus find in the following 
sections that the best fit to both numerics and experiment re- 
quires slight adjustments to the dumbbell model parameters to 
obtain quantitatively optimal fits. 

We also take this opportunity to caution the reader that the 
'microscopic' parameters themselves are subject to change on 
the level of a few percent as experiments and their detailed nu- 
merical modeling evolve (and, hopefully, improve) over time. 
Such changes can be innocuous (e.g. a 1% change to the dia- 
mond lattice constant) but since some of the resulting physics 
is rather delicate, they can feed through to relatively larger 
corrections, most prominently as a factor 3 in the estimated 
value of Td&^ 

From the pyrochlore lattice constant a = 3.54 A one ob- 
tains the diamond lattice constant a^ = ■\/3/2a = 4.34 A. 
Combined with the spin magnetic moment p = lOpB (pb = 
9.27 lO^^"' J/Tesla), this gives the monopole charge q ~ 



As a final remark, it is interesting to compare the force be- 
tween two monopoles at nearest neighbour distance. 
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to that between two eletrons at the same distance, F^i ~ 
1.22 10^^ N, four orders of magnitude stronger! By contrast, 
a pair of Dirac monopoles would experience a force of almost 

IQ-^N. 



B. Internal consistency: screening length vs. monopole 
separation and lattice constant 

The Debye screening length i^Dcbyc is given by the inverse 
of the constant k in Eq. ( 12.101 ). In units of the diamond lattice 
constant ad this amounts to 
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The dependence of Coobye/ad on temperature, after substi- 
tuting p{T) from the numerical solution of Eq. ( 12.131 ) is illus- 
trated in Fig.[3](using for instance A = 4.7 K). 

We anticipate here that there is a systematic discrepancy 
between the DH approximation and the MC simulation results 
on the heat capacity for T > 1 K (see Fig.|7j. To understand 
this, we note the following. 



C. Role of the magnetostatic contribution 
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FIG. 3. Plot of the Debye screening length vs temperature, using the 
density from the numerical solution to the Debye-Hiickel calculation 
in SecHTBl 



Firstly, above T ~ 1 K the screening length becomes 
shorter than the lattice spacing. This artefact arises because 
the DH term in the free energy was derived in the contin- 
uum. For T > 1 K one thus needs to consider the DH results 
with caution. Having said this, once the screening length gets 
very short, the long range nature of the Coulomb interaction 
becomes less important. One can then reliably truncate the 
interactions to short range and use alternative approaches to 
compute the free energy and other thermodynamic quantities, 
as illustrated for instance in AppendixlAl 

Secondly, as T approaches the Curie-Weiss temperature of 
about 2K, the average separation between monopoles, d ^ 
o-d P^^^'^, becomes comparable to the lattice constant a^ and 
the monopole picture is no longer appropriate to describe spin 
ice - monopoles are useful as long as they are sparse, other- 
wise it is more efficient to work directly with the microscopic 
spin degrees of freedom. (In addition, for even higher val- 
ues of T, the neglect of doubly-charged monopoles becomes 
problematic.) For instance, it would be more appropriate to 
use a conventional high-temperature series expansion. 

Another parameter of physical relevance is the ratio of 
screening length to monopole separation: the larger this ra- 
tio, the more appropriate a continuum description is. The di- 
mensionful monopole density pv can be expressed in terms 
of the monopole density per tetrahedron p (which appears 
in the DH calculations in Sec. |II]i using the relation pv = 
3V3p/(8a^). From it, we can obtain the average monopole 

— 1/3 

separation py .By comparing these two length scales, one 
observes that DH theory is near an 'internal' limit of validity, 
as the ratio ■^Dobyo/Py is close to one throughout the range 

of interest. Indeed, ^Dcbyc/Pv ^ ^ ^'^^y below 300 mK, 
dropping by a factor three towards its minimum at 1 K (not 
shown). 



It is interesting to quantify how big the change brought 
about by the DH accounting of Coulomb interactions and 
screening actually is. To do this, let us consider the density 
of monopoles, which will play a role later in the comparison 
with Monte Carlo simulation results (Sec. lIV Al l. In Fig.|4]we 
plot the ratio of the monopole densities from Sec. IIIBI with 
and without the magnetostatic contribution Eq. ( 12. 9t , using 
parameters appropriate for spin ice Dy2Ti207. Within the re- 
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FIG. 4. Ratio of the monopole densities from Sec. lIIBI obtained with 
and without the Debye-Hiickel magnetostatic contribution, Eq. l |2.9t , 
as a function of temperature. 

gion r < 1 K, one notices that DH theory can lead to a more 
than two-times larger monopole density. Given that spin ice 
materials are prone to falling out of equilbrium at tempera- 
tures T < 0.5 K, the behaviour of the system in the temper- 
ature window where DH corrections are sizeable is of crucial 
relevance to experiment. In the limit of low temperatures, the 
DH coiTection instead becomes less and less important. 



D. Monopole-antimonopole pairing 

Debye-Hiickel theory neglects the association of 
monopoles into neutral dipolar pairs (see Ref. [3lj and 
references therein). Although this can in general lead to size- 
able discrepancies between DH predictions and experiments, 
we argue hereafter that pairing corrections are small for the 
observables in spin ice that we consider here, due to the com- 
bination of its limit of validity (T < 1 K, see Sec lIIIBl i and 
the relatively larger energy cost for a monopole excitation, 
A ^ 4 — 5 K, in comparison to the Coulomb energy when 
hard core charges come into "contact" (nearest-neighbour 
distance), £;„„ ~ 3.06 K. 

In order to show this, let us assume that monopoles in spin 
ice are either free (density p^), if separated by a distance larger 
than £b, or bound in a pair, if separated by a distance d shorter 



than £b- Here we choose Is to equal the Bjerrum length, at 
which the thermal energy fc^T equals the Coulomb energy: 



is/ad 



Airaa ^ksT ~ T[K\ 



forDyaTiaOT, (3.6) 



We now consider only Coulomb interactions amongst free 
monopoles and between the two monopoles belonging to the 
same pair, while we neglect monopole-pair and pair-pair in- 
teractions, on the grounds that they are generally weaker and 
they decay faster with distance. We also neglect excluded 
volume effects (therefore, any results we obtain ought to be 
treated with care as the density of monopoles approaches 
unity, which is anyway not the regime we are interested in). 

The free energy /o for the fraction of free monopoles in 
the system is straightforwardly given by Eq. 12.121 The po- 
tential energy term for the bound pairs, of densities pd, d = 
1,2,..., Ib, is also immediate to write as it involves only the 
inter-pair Coulomb term: (2A — Ed)pd, where Ed ^ E^a/d. 
The entropic contribution to the free energy of a bound pair of 
characteristic distance d can be computed from the numbers 
of ways that such pair can appear on the lattice, 
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where Vd is the number of configurations that the two 
monopoles in the pair can take, given say that the centre of 
mass of the pair is fixed. For a nearest-neighbour pair, wi = 2. 
For large values of d, we expect Vd to scale as 2 x And'^. In 
practice, we shall approximate 
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Combining these results, we obtain the free energies (per 
tetrahedron) for free and bound pairs. 



./o = 
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+ T [po In (po/2) + (1 - po) ln(l - po)] (3.10) 
fd = (2A - Ed)pd 

+ T [pd In (pd) + (1 - pd) Hi - pd)] 
-TpdHvd), (3.11) 

as a function of the densities po and pd, d = 1, . . . ,iB- The 
equilibrium free energy of the entire system is then obtained 
minimizing the sum 

./tot = /O + /l + ■ • • + fig 

with respect to po, pi, . . . , pi^ . 

Unlike po, already considered in Sec. IIIBI the pd are ob- 
tained straigthforwardly as 
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Clearly, an intrinsic limit of validity of the theory is given by 
the condition that 



Ptot = Po 
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(3.13) 



In addition, we are of course in particular interested in 

Po > 2 X^dti Pd == Pb- 

The behaviour of po, pi, pb and ptot as a function of tem- 
perature in the regime of interest to spin ice is shown in Fig. |5] 
While at T = 1 K the bound pairs make up for approximately 
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FIG. 5. Behaviour of po (red), pi (blue), pt (green), and ptot 
(open black circles). In the regime of interest to spin ice physics, 
the total monopole density is dominated, at equilibrium, by the free 
monopoles. 

16% of the monopoles in the system, this quickly drops to 7% 
at T ^ 500 mK and to < 10"^% for T < 100 mK. 

Of course, all the considerations in this section apply when 
the system is in thermal equilibrium. This is known not to 
be always the case in experimental settings involving spin ice 
materials! For example, as discussed in Ref. 23, fast varia- 
tions in the temperature of a sample can lead to a "population 
inversion", whereby a relatively high density of monopoles 
survives out of equilibrium down to very low temperatures, 
mostly forming nearest-neighbouring pairs (ptot — Pi)^^- 

The arguments presented in this section are akin to the so- 
called Bjerrum correction to DH. The latter typically leads, 
at low temperatures, to the condensation of all monopoles 
into bound pairs. This is an artifact due to the neglecting of 
monopole-pair interactions, as discussed in Ref.lJU 

Our results do not exhibit any such condensation. The rea- 
son for this difference in behaviour are to be found in the large 
monopole cost with respect to the Coulomb energy at nearest- 
neighbour distance. The net energy gain in the formation a 
bound pair is insufficient to compensate for the corresponding 
entropy loss. The situation would be dramatically different if 
the creation cost of the monopoles were lowered such that it 
can be offset by the Coulomb attraction to another monopole. 



For completeness, we mention that for sufficiently large 
Coulomb attraction the chemical potential of a bound pair 
would have the opposite sign with respect to that of a free 
monopole, leading to a collapse of the system into an ionic 
crystal of monopoles. In spin language, this tranlates into an 
instability of spin ice to an ordered ground state. 



IV. COMPARISON OF DH WITH MONTE CARLO 

We compare the DH results above with Monte Carlo (MC) 
simulations using the spin ice parameters in Ref . |23, reported 
in the previous section. The Ewald summation technique 
was used for the long range dipolar interactions between the 
spins'. We used systems of size IQL^ — 3456 spins (L = 6) 



and single spin flip updates. 



A. Monopole density 

A first comparison between the non-interacting limit and 
the DH approach can be done by looking at the resulting 
monopole density as a function of temperature, Eq. ( 12.51 ) and 
the numerical solution to ( 12.13b . illustrated in Fig. |6] together 
with the monopole density from Monte Carlo simulations of 
dipolar spin ice. 



sensitive to such short range details, such as A, this 8% dis- 
crepancy is not unreasonable. 



B. Heat capacity 

Given the DH free energy (expressed in units of degree 
Kelvin per Dy ion), one can obtain the heat capacity of the 
system in units of J mol^'^K^'^ via the thermodynamic rela- 
tion 
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-NAkBTdl{F/N,kB), 



(4.1) 



where Na is Avogadro's number, j3 = l/ksT, and ks is the 
Boltzmann constant. 

In MC simulations, cy can be obtained by the usual 
fluctuation-dissipation route, measuring the average energy 
(e) and its fluctuations. 
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A comparison between the non-interacting calculations, 
Eq. ( |23]l and Eq. ^^, the DH calculations, Eq. (ITTJl ) 
and Eq. (|2.12t , the single tetrahedron approximation in Ap- 
pendix |Al and Monte Carlo simulations is shown in Fig. |7] 
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FIG. 6. Monopole density from numerical simulations (green tri- 
angles), compared to the analytical result in the non-interacting ap- 
proximation (dashed red line) and in the DH approximation (solid 
blue line). Note that there are no fitting parameters. An improved 
agreement between the simulations and the DH approximation ob- 
tains if we adjust the bare monopole cost to Amc = 4.7 K (black 
dotted curve). 

The agreement between DH and MC results is aheady quite 
reasonable yet it improves considerably if we tune the bare 
monopole cost to Amc = 4.7 K. As mentioned above, we be- 
lieve the origin of this adjustment to be in the short-distance 
physics beyond the dumbbell model of Ref. 12. In quantities 



FIG. 7. Heat capacity from numerical simulations (green triangles), 
compared to the analytical result in the non-interacting approxima- 
tion (dashed red line) and in the DH approximation (solid blue line). 
Note that there are no fitting parameters. Like for the density (cf. 
Fig. |6](, improved agreement between the simulations and the DH 
solution is obtained for a bare monopole cost Amc = 4.7 K (black 
dotted curve). The single-tetrahedron approximation discussed in 
Appendix IaI can only be made to agree with the experimental re- 
sults on a very narrow temperature range, even if we use Jeg as a 
fitting parameter (dash-dotted yellow line). 

Consistently with the monopole density results, a compari- 
son of the heat capacity from DH theory and simulations also 
shows improved agreement using Amc — 4.7 K instead of 
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A = 4.35 K. We shall see in Sec. |Vl]that an 8% larger value 
of A with respect to Eq. ( 13.11 ) is also consistent with the com- 
parison between DH theory and experimental results. 

The results in Fig.|6]and in Fig. [7] clearly show that: (i) a 
theory of point-like Coulomb-interacting charges (in particu- 
lar with the improved value of the bare monopole cost) goes 
a long way into capturing the physics of spin ice, much better 
than conventional approaches based on truncated cluster ex- 
pansions of the free energy of the system; (ii) the long-range 
nature of the interactions is necessary for understanding the 
low-temperature properties of spin ice materials. 



can in principle be tuned straightforwardly, e.g. by decreasing 
D at fixed J^ff, as the magnetic monopole charge is propor- 
tional to D, whereas the scale determining the applicability of 
the monopole picture is set by Joff . 

Indeed, for the nearest-neighbour model with D = 0, where 
there is no magnetic monopole charge, one would be consid- 
ering a Coulomb gas with entropic interactions only. Debye 
screening in such a setting has already been considered in two 
dimensions, for the entropic Coulomb gas encountered in the 
square lattice monomer-dimer modeli2^ 



V. ENTROPIC CHARGE: ROLE OF THE UNDERLYING 
SPINS 

In disregarding the underlying spins in the Debye-Hiickel 
approximation to the free energy of spin ice, we fail to ac- 
count for quadrupolar corrections to the monopole descrip- 
tioni^ (of which we have seen an effect in the value of the 
bare monopole cost A). We also neglect additional spin en- 
tropic contributions (other than the entropy of mixing of the 
monopoles)^*^. 

The latter take the form of an entropic charge that adds 
onto the real magnetic charge (or, rather, magnetic and en- 
tropic coupling constants add) for the monopole Coulomb in- 
teractions. In Appendix iBl we derive an analytical expression 
for the entropic interaction strength and confirm the result by 
comparing it to Monte Carlo simulations. One can then repeat 
the DH calculations including the entropic correction. The re- 
sults are shown in Fig. [8] (dashed cyan lines), in comparison 
to the previous results (solid blue lines), for the parameters in 
Sec. |IV]with Amc = 4.7 K. The behaviour of the monopole 



VI. EXPERIMENT 

We now proceed to compare the DH results with experi- 
mental data on Dy2Ti207. We find good agreement, which is 
further improved if we use the latest material parameters from 
Ref. 26 instead of those in Ref.|23. Namely, the magnetic mo- 
ment of the rare earth ions is 9.87 ^b instead of 10 iib\ the 
diamond lattice constant is 4.38 A instead of 4.34 A; and the 
nearest-neighbour exchange coupling varies between —3.53 
and -3.26, instead of J = -3.72 K. 

These values result in a new magnetic monopole charge 
of 4.5 /is/A; a nearest-neighbour interaction strength be- 
tween monopoles £"„„ = 2.88 K instead of 3.06 K; a dipo- 
lar coupling constant D = 1.32 K instead of 1.41 K; and 
a bare monopole cost in the range (4.05,4.23) K instead of 
A = 4.35 K. We reiterate that there are also small correc- 
tions due to further-range superexchange and the quadrupo- 
lar interactions, which are not easily incorporate into the DH 
framework. 



A. Heat capacity 
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FIG. 8. Effects of the entropic charge (dashed cyan lines) on the 
Debye-Hiickel estimate of the heat capacity and monopole density 
(solid blue lines). 

density and of the heat capacity clearly show that the entropic 
contribution can be safely neglected in the low temperature 
regime where the DH approximation is valid. It is worth not- 
ing that the relative strength of magnetic and entropic charges 



A comparison between the experimentally measured heat 
capacity and the one obtained from DH theory, shows again 
that the bare monopole cost A e (4.05, 4.23) K from Eq. ( 1311 ) 
is somewhat too small. Better agreement can be obtained if, as 
in the comparison with MC simulations, we allow for an 8% 
increase in the value of A g (4.37, 4.57) K (see Fig.|9l). This 
is in agreement with the results presented in Ref. [lOi (Fig. 1), 
where a value of A = 4.35 Ki^ was used. 



B. 'Dressed' monopole energy and AC susceptibility 

The bare monopole cost A is half the energy required for 
creating and separating to infinity a pair of monopoles against 
their long-range Coulomb attraction. When other monopoles 
are present, screening effectively truncates the range of the 
interactions and there is no further energy cost to separating 
a pair beyond the screening length. In this case it is more 
appropriate to consider the 'dressed' monopole energy A^ as 
the energy per monopole that it takes to create a pair and sep- 
arate it beyond the screening length. It is indeed the energy 
Ac; - rather than A - that controls for instance the equilib- 
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FIG. 9. Experimental results for the heat capacity of Dy2Ti207 
(black squares) from Ref. [lOl in units of J/mol K, compared to the 
analytical result from Debye-Huckel theory with A — 4.37 K (solid 
blue line) and A = 4.57 K (dashed cyan line). 



rium density of the monopoles p ^ e^^''^'^ at intermediate 
temperatures. 

Given the creation energy for a nearest neighbour pair 
As = 2A — Enn and the expression for the DH screening 
length, Eq. ( 13.51 ). one obtains 



2Ad(T) = 2A - S„ 
= 2A - £:„ 
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(6.1) 



whose behaviour is illustrated in the inset of Fig. [TO] 

A place where this screening effect of the magnetic 
monopoles becomes particularly evident is in susceptibility 
measurements of magnetic relaxation time scalesi^^^l. Given 
that the monopoles are responsible for any changes in mag- 
netisation in a spin ice configuration, the ability of the sys- 
tem to respond to an applied magnetic field is affected by the 
monopole density. For non-interacting monopoles, Ryzhkin 
showed that in the low temperature, hydrodynamic regime the 
characteristic susceptibility time scale r is inversely propor- 
tional to the monopole density^S, 



r ^ oc i/Tp{T), 



(6.2) 



where i/ is the mobility of the monopoles. This result is Ukely 
to be asymptotically correct as T — ?► at zero wavevector even 
in presence of Coulomb interactions, although it is modified 
at finite wave vectors. 

In App. ICl we show that i' ^ 1/T under the assumption 
that Metropolis dynamics are a good approximation to the mi- 
croscopic spin flip processes in spin ice. Therefore, 



Tcxi/p(r). 



(6.3) 



As we argued above, at intermediate temperatures p{T) is 
controlled by the dressed monopole energy Ad{T) rather than 



FIG. 10. Experimental magnetic relaxation time scale r as a function 
of temperature from susceptibility data, Ref.[l9l (black open squares). 
The rapid increase in r at low temperatures is due to the paucity 
of defects responsible for the magnetic rearrangement of a spin ice 
configuration (namely, the monopoles). This increase cannot be de- 
scribed by a single exponential (activated behaviour), as it is evi- 
dent for instance by comparison with the curve r = ro exp(A/T) 
(dashed magenta line), say with A — 4.5 K. On the contrary, a 
much better agreement is obtained if we replace the bare monopole 
energy A with the 'dressed' energy Ad(T) (solid blue curve for 
A = 4.37 K and solid cyan curve for A — 4.57 K). This is 
compared to r ex 1/p, where p is obtained from the DH approx- 
imation (blue open circles for A = 4.37 K and cyan open circles 
for A — 4.57 K), showing that indeed the dressing of A accounts 
for the leading non-exponential correction in the temperature depen- 
dence in the monopole density. The microscopic time scale was set 
by imposing that the analytical results pass through the experimen- 
tal data point at 4 K (see Ref. 21). The inset shows the 'dressed' 
monopole energy Ad as a function of temperature (solid blue curve 
for A — 4.37 K and solid cyan curve for A = 4.57 K). 



the bare energy A. Indeed, r is poorly fitted by a single ex- 
ponentiali^i^ such as r = tq exp(A/T). On the contrary, the 
curve T — Tq exp[Arf(r)/r], captures correctly the faster- 
than-exponential grows of t at low temperatures, despite the 
fact that it still significantly underestimates the experimental 
value of T (see Fig.fTOb.^" 

Given the good agreement between DH theory and experi- 
ments regarding the heat capacity of the system (Fig. |9]) and 
given that a similarly good agreement in the heat capacity 
from Monte Carlo simulations implied a good agreement also 
for the monopole density (Fig. |6] and Fig. |7]), one would ex- 
pect that p{T) from Debye-Hiickel used in Fig.[TO]is in fact a 
good estimate of the experimental monopole density. There- 
fore, the fact that Eq. (16.31 ) underestimates the experimental 
results even when using p{T) from DH theory is likely due 
to corrections to the dependence r ex l/p{T) arising from 
Coulomb interactions at intermediate monopole densities. 

At the lowest temperatures (provided of course no ordering 
or freezing intervenes, as it likely would), when monopole 
separation and screening length both diverge, the effective 
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Ad -^ A, and hence we expect the superexponential be- 
haviour to go away and the curve to follow the standard Ar- 
rhenius behaviour r ^ exp(A/r). 

From a purely phenomenological perspective, it is interest- 
ing to notice that a very good agreement beween DH theory 
and experiments on the susceptibility time scale r (at interme- 
diate temperatures) can be obtained by substituting Eq. ( I6.3l l 
with T ex l/p''(r), with 7] = 3/2 for A = 4.37 K and 
7] — 4/3 for A = 4.57 K (see Fig. [TTT i. Further work is 
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FIG. 1 1 . Experimental magnetic relaxation time scale r as a func- 
tion of temperature from susceptibility data, Ref. 4? (black open 
squares). The temperature dependence is captures very accurately 
by a phenomenological equation of the type r oc 1/p^, where p 
is obtained from the DH approximation (blue upward triangles for 
A = 4.37 K and rj — 3/2; cyan downward triangles for A = 4.57 K 
and rj — 4/3). The dashed magenta line illustrates the curve 
T — To exp(A/r) with A = 4.5 K, for comparison. 

needed to understand the reasons behind such a good overlap. 



in the phase diagram of the system-", with an exponential in- 
crease in the low-temperature fraction of neutral pairs drain- 
ing the free monopole density to zero. This can (and ought to) 
be compensated by a further extension to include interactions 
between dipolar bound pairs and free monopoles, leading to 
the so called dipole-ionic (DI) contribution-' k The full DHB- 
jDI theory indeed cures the unphysical features identified for 
DHBj, while remaining of course only an approximation to 
the exact free energy of the system. 

Further improvements on the DHBjDI theory include ac- 
counting for hard-core (HC) effects^. It is certainly worth- 
while developing the theory further in this direction, espe- 
cially in settings or for quantitites where new phenomena 
(e.g., a dominant population of bound pairs), rather than only 
quantitative corrections, ensue. 



VIII. CONCLUSIONS 

In summary, we have presented a theory for the low- 
temperature physics of spin ice within the Debye-Hiickel 
framework familiar from the study of (electric) Coulomb liq- 
uids. The success of this simple approach in treating the low- 
energy physics of spin ice is a testament to the power of the 
'variable transformation' from magnetic dipoles to magnetic 
monopoles appropriate to the Coulomb phase with its emer- 
gent gauge field. 

With this first step accomplished, next on the wishlist are 
a number of items some of which should push our atten- 
tion beyond the framework provdided by the DH paradigm. 
Firstly, a more detailed understanding of spin ice (hydro- 
)dynamics; secondly, an extension of this theory to a broader 
class of parent Hamiltonians, perhaps even including coher- 
ent quantum dynamics; and thirdly, contact with all the non- 
equilibrium experiments suggesting that not only the sparse- 
ness of monopoles but also phononic physics plays a role in 
the freezing of spin ice around Tf^ 



VII. BEYOND DEBYE-HUCKEL 

Debye-Hiickel theory is probably the simplest approxima- 
tion to obtain the free energy of a gas of Coulomb interacting 
particles short of ignoring interactions altogether. A number 
of improvements are available in the vast literature on the sub- 
ject"", which one can use to obtain a more accurate description 
of the magnetic monopole behaviour in spin ice. 

Without actually implementing them, we briefly recall 
hereafter two common extensions of the DH model. Firstly, 
Debye-Hiickel theory neglects the association of monopoles 
into neutral dipolar pairs, which we have already briefly dis- 
cussed above (see Ref. [31] and references therein). Following 
Bjerrum^^ (Bj) one can account for such bound pairs, thus 
compensating in good part for the uncontrolled linearisation 
of the Poisson-Boltzmann equation that is at the basis of the 
DH self-consistent solution. However, whilst being an overall 
refinement of DH, DHBj theory leads to unrealistic features 
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Appendix A: Single tetrahedron approximation 

An alternative approximation that can be used to obtain the 
spin ice free energy and related thermodynamic quantities is to 
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use a truncated cluster expansion. Most simply, this amounts 
to computing explicitly the free energy of an isolated tetrahe- 
dron by direct summation over all 2^ states. 

At this level, all interactions are nearest-neighbour ones. In 
terms of this effective short range coupling Jeff, the partition 
function of a tetrahedron is 



6 + Se"^'^''"/^ + 2e"^''"ff/^ 



Nt 



(Al) 



From this, one can estimate the partition function of the entire 

system. 



Z = 2 



N, 



6 + 8e-2^<=«/^ + 2e-8''=«-/^ 
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, (A2) 



and thus the free energy per spin in degrees Kelvin, 

F/N,kB = -iT/N,)lnZ. 

Substituting into Eq. (14.1b . we obtain the heat capacity of 
the system (in units of J/K per Dy ion). 
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The choice of Jeff — 5D/3 + J/3 = 1.11 K, which corre- 
sponds to the nearest-neighbour interaction strength from the 
exchange plus dipolar coupling constants, yields a very poor 
agreement with the experimental data (not shown). The sit- 
uation improves slightly if we take advantage of the projec- 
tive equivalence between dipolar and nearest-neighbour in- 
teractions on the pyrochlore lattice-'"*. Instead of truncating 
the dipolar contribution to 5D/3, one can therefore use the 
effective value of Jn„ that yields the same low-energy spec- 
trum as from the long range dipolar interactions. This value 
can be derived using the dumbell decomposition in Ref. Jjl 
Jeff = 1.45 K. The result is shown in Fig. 1 of Ref. llOl and 
it is indeed in quantitative agreement with the experimental 
data at high temperatures T > 2 K, as expected of a cluster 
expansion of the free energy. 

Note that even if we allow Jeff to vary as a fitting parameter 
in the theory, the shape of cy {T) does not change significantly 
and it can be brought to agree with the experimental data only 
over a very narrow temperature interval. By comparison, this 
highlights even more how effective the Debye-Hiickel free en- 
ergy is at capturing the low energy fluctuations in dipolar spin 
ice. 



Appendix B: Entropic monopole charge 

The effective description of spin-ice in the absence of 
monopoles is given by the probability distribution of a 
magnetostatic-like (divergenceless) field^^ 
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The first term Eq. ( IBll i is purely entropic in origin. The ge- 
ometric field i?°"' (r) is obtained from coarse graining fixed- 
length vectors that identify the local direction of the spins in 
the system. Here Vcci\ is the volume of the primitive unit cell 
(Fig.[T2]i. Introducing the coarse grained (dimensionless) field 
B{r) defined at the centre of each tetrahedron (belonging to 
one of the two sublattices) as 
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the stiffness coefficient can be determined to be /C = 3/8. 
(Note that we used a different field normalisation with respect 
to Ref. ,351 so as to preserve the underlying spin length equal 
tol.) 




^"l = 4/xf3a. 



FIG. 12. Lattice conventions. The highlighted portion of the blue 
cube (i.e., the 16-spin cubic unit cell in spin ice) corresponds to a 
possible choice of the primitive unit cell in the fee lattice formed by 
the centres of one sublattice of tetrahedra in the pyrochlore lattice 
(circled in green in the figure). 



The second term Eq. (lB2b accounts for the magnetic energy 
stored in a spin ice configuration (devoid of monopoles). In 
this case, H'^^^^{r) is the magnetic field generated by the spin 
magnetic moments /i pointing in the local spin direction (/io is 
the permeability of the vacuum, ks is the Boltzmann constant 
and T is the temperature of the system). 

Given that the total field B = ijlq^H + M) is always di- 
vergenceless, the field H^^^{r) can be equivalently replaced 
by the magnetisation per unit volume M, which in turn can 
be obtained by coarse graining the spin magnetic moments. 
Using the scheme (IB3b already adopted for B°'^^{r) over a 
primitive unit cell, we have that 



7J'"^«(r) = M{r) 
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Therefore, the difference between the two terms Eq. ( IB II ) and 
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Eq. ( IB2I ) can be reduced to different coefficients 
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to the same integral / |B''"*(r)p d^r. 

It is convenient to re-express the magnetic coefficient in 
terms of the magnetic Coulomb energy of two monopoles 
placed in adjacent tetrahedra (expressed in degrees Kelvin), 
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where we used the fact that q = 2/i/ad, a^ being the diamond 
lattice constant. By comparison with the entropic coefficient, 
we can then identify the entropic counterpart to the neareast- 
neighbour Coulomb energy. 
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If we finally use the fact that Ucoii is 1/4 of the volume of 
the 16-spin cubic unit cell in spin ice, v ~ (Aad/V^)^, and 
that with the coarse graining ( |B3l l /C = 3/8, we arrive at the 
result 
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It is interesting to convert this value into an entropic 
monopole charge : 
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The entropic charge of a monopole becomes larger than the 
real magnetic charge only for T > 8 K, well beyond the limit 
of validity of the monopole description of spin ice. In the ex- 
perimentally relevant temperature range 0. 1 — 1 K, the entropic 
contribution ranges from 1% to 10% of the real magnetic con- 
tribution to the energy of the monopoles. 

In order to confirm this analytical estimate of the entropic 
Coulomb interaction strength in spin ice, we have run Monte 
Carlo simulations of the nearest-neighbour spin ice model, 
sampling only configurations with two monopoles (one pos- 
itive, one negative). Such configurations are all isoenergetic 
and the monopole positions can be updated at every Monte 
Carlo step without rejection. Ergodicity was tested by com- 
puting spin-spin autocorrelation functions. The distribution 
of separation distances between the two monopoles was then 
sampled both in Monte Carlo time and across different initial 
configurations and random number seeds. 

From Eq. IB II it follows that the entropic interaction be- 
tween the two monopoles leads to a probability distribution 
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FIG. 13. Top Panel: Distribution of distances per lattice site be- 
tween two monopoles in a spin ice configuration of 16 x L* spins, 
L = 64. (top panel, red curve). The expected form due to 
the entropic Coulombic interaction is P ~ exp{E'^'^ /TR) and 
the solid yellow line is the linear fit of In P{R) as a function of 
1/R. Bottom Panel: Finite size scaling of the nearest neighbour 
entropic interaction E"^ /T vs. the inverse system size 1/L, L = 
16, 32, 48, 64, 80, 100. The dashed black line and shaded cyan re- 
gion are a guide to the eye for a reasonable L — >■ oo extrapolation 
and confidence interval, leading to El'^ /T ~ 0.375 ± 0.015. 



of the form V{R) - R^ exp(£:™7ri?), where R is the sep- 
aration distance in units of the diamond lattice spacing. In 
particular, if we sample the distribution per lattice site at dis- 
tance R, it has a purely exponential form ^ exp(_E™*/Ti?), 
and one can obtain the value of E°^^ /T from linear fits in 
semi-logarithmic scale (Fig. [13] top panel). 

We repeated these fits for different system sizes in order 
to account for finite size scaling (illustrated in Fig. [T3j bot- 
tom panel). Even though the accuracy of our simulations 
does not allow for a reliable extrapolation in the L — > oo 
limit, the nearest-neighbour entropic interaction strength ap- 
pears to lie in the interval £^™V^ - 0.375 ± 0.015, in rea- 
sonable agreement with the analytical value in Eq. JBIOI ), 
2/V^7r ~ 0.36755. 
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Appendix C: Monopole mobility 

The mobility of the monopoles in spin ice (and thus its tem- 
perature dependence) can be estimated from microscopic con- 
siderations, under the assumption that Metropolis-like equa- 
tions govern the dynamics of the system^-'-'*. 

The mobility of a particle is given by the ratio of its drift 
velocity Vd over the driving force strength qE, v = Vd/{qE). 

Under Metropolis dynamics for a particle with charge q in 
a field E, the average displacement in a single step is 



Aa;^^ 
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where i is the characteristic microscopic length scale, V is the 
potential difference for a single hopping process, 1 is the prob- 
ability to hop in the direction of the field, and e'xp{—j3qV) is 
the probability to hop in the opposite direction. 

Note that, on a lattice, there can be several inequivalent for- 
ward and backward hoppings, depending on the direction of 
the field. For example, while a 45° field applied to charged 
particles living on a square lattice is described straightfor- 
wardly by the above equation (with C — a/\/2, a being the 
lattice spacing), a 90° field on the same lattice allows for a for- 
ward, a backward, and two perpendicular hopping processes 
(see Fig.O. One therefore needs to average over all of them 





FIG. 14. Two examples of how the available hopping processes 
depend on the direction of the applied field on a square lattice: a 45° 
field (left) and a 90° field (right). 



to obtain the correct value of Ax. 

For convenience, we choose to define the mobility i' as 



Ax/a 



e 1 



-t3qV 



To 



flTo 1 + e-l^i^ 
V qEa, 



(C2) 
(C3) 



for small values of the applied field E. Here a is the (dimen- 
sionful) lattice constant and tq is the microscopic time scale 
for a single MC step. At large temperatures with respect to 
the field strength, one can expand the exponentials and arrive 
at the expression 



1 1 



-HqV 



To qEa al + e-'^?^ 
1 eV/{Ea) 



Tq a 2kBT 



O 



qEa 



(C4) 
(C5) 



For example, the case of a generic field direction on the 
anisotropic square lattice, with lattice constants a and b, gives 



1 1 acosO + hsine -aco&ee~^'i^'"'°'''' -hsuvOe-P'i^^' 



To qEa"^ 



To 'iksT 



1 



62 -"2 



1 

sin' 



f^ — j3qEa cos I f,—l3qEbsinf: 

- + {/3'E) . 



(C6) 



If the lattice is isotropic (a — b) the mobility is independent 
of the direction of the applied field. 



1 1 



O {/3^E) 



(C7) 



To iksT 
The mobility of monopoles on an isotropic diamond lattice. 



I 

of lattice constant ad, with respect to a generic field direction 
e can be computed in a similar way, with the additional care 
that there are now two inequivalent sublattices. With respect 
to one sublatice, we obtain 
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I I 1 (ei + 62 + es) min 


I g/3qEad(ei+e2+e3)/\/3 


+ (ei- 


- 62 - 63) min 


1 g/3gEa 


diei-e2-e3)/-/3 


To qEaa V3 min 


l^ gl3qEad{ei+e2+e3)/V^ 


+ min 


l^gl3qEaa{ei-e2 


-&3)/V3 





+ (-61 + 62 - 63) min 



2 gPqEad(-ei+e2-e3)/\/3 



2 g/3(ir-Ead(-ei+e2-e3)/y3 



+ (-61 - 62 + 63) min 



2 g/3g-Eajj(-ei-e2+e3)/y3 



2 p0qEad{~-ei-e2+e3)/\/3 



[( 



61 



62 + 63)^ 9<(ei + 62 + 63) + (ei - 62 - 63)^ 0<(ei - 62 - 63) 



To 12kBT 

+ (-61 + 62 - 63)2 e<(-ei + 62 - 63) + (-61 - 62 + 63)^ e<(-ei - 62 + 63)] + O {/3^ E) 

where 0< (a;) — 8 (—a;) is the Heaviside theta function. With respect to the other sublatice, we obtain 

[(ei + 62 + 63)^ [1 - 6<(ei + 62 + 63)] + (ei - 62 - 63)^ [1 - 6<(ei - 62 - 63)] 



(C8) 



To UkeT 

+ (-61 + 62 - 63)^ [1 - e<(-ei + 62 - 63)] + (-61 - 62 + 63)^ [1 - e<(-ei - 62 + 63)]] + O {(5^ E) 



(C9) 



If we finally take the average of both sublattices, we arrive at 
1 1 



2to l2kBT 



[(ei + 62 + 63)^ + (ei - 62 - 63)'^ 



{P'^E) 

1 1 1 



6 To ksT 



+ (-61 +62 - 63)^ + (-61 



O {l3^E) 



62 + 63J 



(CIO) 



independently of the direction of the field E. 

If the magnetic monopoles on the diamond lattice are in 
fact the collective excitations in a spin ice system, one needs 
to take into account the constraint that one of the three pos- 
sible hopping directions is essentially forbidden, as it would 
create doubly charged excitations. Taking the average over the 
possible forbidden directions does not introduce a dependence 



I 

on the field direction and we can therefore choose to compute 
the mobility in a [100] magnetic field for convenience: 



1 1 



1 2 



-|3qE^aa/^/3) 



2to qEaa ^3 2 + e-l^iEiad/V3) 
1 1 1 1 - 26"'^''-^("''/^) 

2^ qEud V3 1 + 2e-/59-E(«<i/^) 
4 1 1 



27 TnkBT' 



(Cll) 



Notice that these results are independent of whether the po- 
tential and field had an entropic or magnetic origin, provided 
that the assumption of the field being smooth over distances 
of the order of the lattice spacing a^ holds. This definition of 
the mobility shows in fact that it depends only on some mi- 
croscopic time scale tq and on the thermal energy per particle 
in the system. 
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